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Abstract. We present the results of a numerical simulation of propagation of cosmic 
rays with energy above 10^^ eV in a complex magnetic field, made in general of a large 
scale component and a turbulent component. Several configurations are investigated 
that may represent specific aspects of a realistic magnetic field of the Galaxy, though 
the main purpose of this investigation is not to achieve a realistic description of the 
propagation in the Galaxy, but rather to assess the role of several effects that define 
the complex problem of propagation. Our simulations of Cosmic Rays in the Galaxy 
will be presented in Paper II. We identified several effects that are difficult to interpret 
in a purely diffusive approach and that play a crucial role in the propagation of cosmic 
rays in the complex magnetic field of the Galaxy. We discuss at length the problem 
of the extrapolation of our results to much lower energies where data are available 
on the confinement time of cosmic rays in the Galaxy. The confinement time and 
its dependence on particles' rigidity are crucial ingredients for 1) relating the source 
spectrum to the observed cosmic ray spectrum; 2) quantifying the production of light 
elements by spallation; 3) predicting the anisotropy as a function of energy. 
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1. Introduction 

A complete understanding of tlie origin of Cosmic Rays (CRs) will be achieved when the 
acceleration processes, the sources and the propagation from the sources to the Earth 
will be included in a self-consistent theoretical framework. This goal is far from being 
achieved: for ultra high energy cosmic rays (UHECRs) the issue of the propagation is 
probably easier to understand since the effect of extragalactic magnetic field is expected 
to be not crucial, at least above energies of ~ 4 x 10^^ eV [Ij. On the other hand in this 
case the sources are fully unknown. For CRs below ~ 10'^ — 10^ GeV we are confident 
that the sources are located within our Galaxy and most likely are supernova remnants 
(SNRs) [2]. In this respect a large bulk of information is being collected from X-ray 
and 7-ray astronomy: high resolution X-ray observations have shown the presence of 
intense magnetic fields in the vicinity of the shocks that bound the shell of the remnant 
[3], thereby making the acceleration process easier to understand. The combination 
with the observed X-ray spectra and the outstanding detection of 10 — 100 TeV gamma 
rays from a few SNRs [1] make a rather strong case in favor of these astrophysical 
sources being the accelerators of protons up to the knee or slightly above it [2]. Nuclei 
would then be accelerated to even higher maximum energies, up to ~ 10^^ eV for 
iron nuclei. Although the observational situation and the theoretical understanding are 
both experiencing a substantial improvement as far as the sources (or at least SNRs) are 
concerned, a realistic description of the propagation of CRs in the interstellar medium 
(ISM) is still missing, despite the very impressive amount of work carried out on the 
topic (see [5] and references therein for a recent review). Such work may be classified 
in two broad classes: analytical approaches and simulations. 

Most analytical work is based on the solution of the diffusion-convection equation 
from a distribution of sources in a medium with given diffusion properties. We include in 
this class the work that is based on a numerical solution of the transport equation (e.g. 
GALPROP [6] or the model presented in [7]). In the most general case, the equation has 
been solved with both parallel and perpendicular diffusion taken into account. These 
approaches start from the premise that the magnetic field of the Galaxy induces only 
a diffusive behaviour on CRs, namely the turbulent field is the key ingredient. This 
component is provided a priori, either in the form of pre-calculated diffusion coefficients 
or in the form of turbulent spectra. It is worth stressing that the spectrum of the 
turbulence responsible for particle diffusion, the total power in turbulent modes and the 
origin of such turbulence are unknown. However, if one assumes that the spectrum is 
known, then the diffusion coefficients could be calculated, at least in principle, using 
quasi-linear theory and neglecting the geometry of the large scale background magnetic 
field. 

The most common approach in the literature consists of using low energy data on 
the ratio of secondary to primary nuclei in CRs as a function of energy to infer the 
energy dependence of the propagation time, which in turn leads to a rough knowledge 
of the energy dependence of the diffusion coefficient. Such dependence is then adopted 
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in the solution of the transport equation. 

The shortcomings and advantages of using the diffusion equation to describe the 
propagation of CRs are easy to identify: this approach allows one to achieve a basic 
understanding of some issues (for instance the spectral steepening induced by the particle 
propagation and leakage from the Galaxy). Moreover the approach can be used without 
limitations in the dynamical range (particles from very low to very high momenta can 
be included). On the other hand, the diffusion coefficients are given quantities; even 
when the diffusion coefficients (as functions of particle momentum and spatial location) 
are calculated from first principles (quasi-linear theory) they are often used in regimes 
where the initial assumptions do not necessarily hold. In addition, a multitude of effects 
related to spatial gradients of the large scale fields are hardly accounted for. 

The numerical simulation of the propagation of CRs in arbitrary magnetic fields 
solves part of these problems, but is limited by the constraints on the computational 
time. Previous investigations using this technique concentrated on very high energy 
cosmic rays (~ 10^^~^^eV), and on the deflections produced by their passage into the 
Galactic Magnetic Field [8] or on their anisotropy around 10^^ eV |9]. Other attempts 
investigated lower energies, e.g. Ref. [10] was able to reach down to 10^'' eV, and 
calculated the times of escape from the galaxy as a function of energy. The results 
obtained, however, seem to be inconsistent with measurements at low energy. Indeed, 
in Ref. ^Q\, the escape time at 10^^ eV is found to be of the order of 10^ yr with an 
energy dependence of E~^, much steeper than the one expected for example from a 
normal Kolmogorov turbulence. The extrapolation of this value to 10^ eV produces a 
value several orders of magnitude larger than the measured one. The problem of the 
steepness of the escape time in the simulations seems to be a general one: it is present 
also in our simulations and seems to continue to lower energies as we discuss below. 

In this paper we describe the numerical code that we recently completed for the 
propagation of cosmic rays in arbitrary magnetic fields (both in their large scale and 
turbulent components). The code represents a substantial improvement on previous 
efforts in the same direction in several ways: first, we succeeded in propagating the 
particles down to energies of 10^^ eV, lower by at least one/two orders of magnitude 
compared with previous simulations. Second, the turbulence responsible for diffusive 
particle motion can be taken as three dimensional or one-dimensional, and as isotropic 
or anisotropic. The large scale field is also completely arbitrary. 

We present the results of this simulation effort in two papers. In the present paper 
(Paper I) we discuss all technical aspects and apply the approach to several toy models 
of the magnetic field of the Galaxy in order to emphasize the role of the physical effects 
that is necessary to include in order to understand the propagation of CRs. In a second 
paper (Paper II) we will describe the results of the simulation for given configurations 
of Galactic magnetic fields which are commonly assumed as realistic. 

The paper is organized as follows: In ^ we describe the technical aspects of the 
simulation, with special attention to the generation of the turbulent magnetic field. In 
^ we illustrate some basic concepts of diffusion in the context of quasi-linear theory. 
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which allows us to define what is the common lore of cosmic ray propagation in the 
Galaxy in terms of diffusion and drifts. In §1] we describe the numerical procedure 
adopted to calculate the parallel and perpendicular diffusion coefficients. Finally, in ^ 
we describe the results of our computations for several toy models of the large scale field 
of the Galaxy. We present our conclusions in ^ 

2. Description of the simulation 

We propagate particles in a magnetic field, B = 5B + Bq, that is the sum in each 
point of a regular and a random component. Both of them can in principle depend 
on the position. As detailed in Ref. [H] there are basically two methods to implement 
the turbulent field: 1) pre-computing the field on a grid using Fast Fourier Transform 
(FFT) and 2) calculating the field in each point along the particle trajectory as the 
superposition of plane waves [T2] . 

In the FFT approach the field is pre-computed on a grid in real space from its 
power spectrum in reciprocal space. We set up a three dimensional grid with integer 
coordinates from to — 1. Each vertex on the grid corresponds to a wave vector k 
with components given by the coordinates of the vertex. If any one of the components 
of fc, for example kx-, is larger than N/2, then we substitute it with — (A^ — k^) in order 
to take into account negative frequencies. For each k we construct an amplitude vector, 
-Bfc, with a length proportional to the square root of the power in the corresponding 
mode: /c"*^"^"*"^)/^, a random direction in the plane orthogonal to k and a random phase. 
Choosing the amplitude proportional to /i;~(T+2)/2 jjiakes sure that the power spectrum of 
the turbulent field is V{k) oc k~^ , whereas choosing the direction in the plane orthogonal 
to k assures that V ■ 6B = 0. We also have to make sure that B^ satisfies the following 
condition for the resulting magnetic field to be real: B(^ki,k2,k3) ~ ^*[N-ki,N-k2,N-kz)- 
The normalization is obtained by requiring that = S-^fc and -Bfc=(o,o,o) is set to 

to have {5B) = 0. At this point we compute the FFT [13] and obtain the turbulent 
field defined on a cubic grid with side Lmax and spacing Lmin = L^^^^/N . We typically 
use A^ = 256. 

We assume the box is replicated periodically all over the simulation volume and 
in order to calculate the turbulent field in a given point the code uses the field value 
of the closest vertex. Another possibility is to do an interpolation of the values at the 
eight vertexes surrounding the point. We verified that the results obtained with the two 
methods are equal on scales larger than the cell size {L^in) and we decided to use the 
former method. 

The above description is valid for the general case of isotropic turbulence. We also 
used ID turbulence, a superposition of Alfven waves, and in this case the generation 
proceeds along the same lines, but the fcs are now chosen only parallel to the background 
field, so that the fiuctuating magnetic field is always perpendicular to it. For the ID 
field we typically use A^ = 4096. 

In the second approach the field is constructed as the sum of A^ plane waves 



Numerical propagation of high energy cosmic rays in the Galaxy I: technical issues 5 



SB = Y^ Ak„en exp{iknz'^ + , (1) 

n=l 

where e„ = cosa^i^ + isina^^^ and a„ and are random phases. The primed 
coordinates are obtained by rotating the reference frame so that the z axis coincides 
with the direction of propagation of the n-th wave, The directions of the Nm waves 
are chosen randomly, while their amplitudes, A^.^, are chosen as a function of |fc„| 
according to the type of turbulence wanted. We follow Ref. [12] and we use: 

Al = a'G{k)[j:Gik^r\ (2) 

n=l 

where 

In these equations a fixes the normalization of the field, o"^ = {SB"^), Lc is the correlation 
length, 7 is the slope of the turbulence power spectrum, d is its dimensionality and AV^'^'^ 
is the volume element for the chosen dimensionality. In the present paper we use 3D and 
ID turbulence and in these cases AV^^^ = Ank'^Ak and AV^^^ = Ak. The wavenumbers 
are chosen evenly spaced in logarithmic scale between kmm and /cmax and Ak = kA log k. 

The number of waves, N^, used in the summation ([1]), is a key parameter and it 
should be large enough to reasonably describe the turbulence. In Ref. [13] it was shown 
that if Nm is too small the transition from rectilinear to diffusive propagation occurs 
on timescales much larger than the correct ones. It was also found that a value of 100 
waves per decade is a reasonable compromise between accuracy and computation time 
and we use this value in our simulations. 

The only difference for the case of ID turbulence is that instead of choosing the fc„ 
isotropically we choose them in the direction parallel to the background field. 

Both the methods described have their advantages and disadvantages: with the 
FFT approach the time needed to obtain the turbulent field in a given point is in 
general much smaller than in the plane wave approach. In the first case all that is 
required is a lookup in a table (and possibly some interpolations), whereas in the latter 
case there is a summation over hundreds of waves to be done. On the other hand, the 
dynamic range of the turbulence, L^^x/ Lmm, in the FFT approach is limited (at least 
in the isotropic turbulence case) by the memory available to store the huge matrices 
describing the magnetic field grid, whereas in the plane wave approach the memory 
limitations are absent and the dynamic range can be as big as required with the only 
limit given by the computation time. As mentioned in Ref. [11] , other limitations of the 
FFT approach are inherent in its discreteness, L^in, and in its limited size, Lmax, and 
the results obtained with it can not be trusted when the Larmor radius of the particles 
is smaller than Lmin or larger than Lmax- 
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3. Basic facts about diffusion and drifts 

In this section we summarize the basic facts on diffusion of cosmic rays in a turbulent 
magnetic field superimposed to a large scale spatially constant magnetic field Bq = BqZ. 
Gradients in the large scale field induce drift motions of the particles that add to the 
diffusive motion and in fact in some circumstances may even dominate upon diffusion. 



3.1. Diffusion 

In all the cases that we consider below we investigate 3D turbulence, namely the 
perturbation of the large scale field has components both in the plane perpendicular 
to Bq and along Bq- Therefore this case is somewhat more complex but supposedly 
more realistic than the simpler case of Alfven waves propagating along the field Bq 
(we refer to this the ID case), typically considered in the literature on quasi- 

linear theory. In the case of 3D turbulence, the perpendicular diffusion, though small 
compared with the parallel diffusion in the quasi-linear regime, may become important 
for the cases of strong turbulence SB/Bq > 1. On the basis of quasi-linear theory the 
ratio of perpendicular to parallel diffusion coefficient is given by 

£i = I (4) 

where the parallel pathlength is Ay = 3D\\/v, v is the particle velocity and is the 
Larmor radius in the unperturbed magnetic field Bq. This expression remains valid as 
long as 6B/Bq <^ 1, but it also suggests that the perpendicular diffusion coefficient 
approaches the parallel diffusion coefficient in the regime of strong turbulence. In fact 
the real ratio of the diffusion coefficients is affected by the random walk of the field 
lines, which is not taken into account in Eq. HI This fact was found in [15] and further 
discussed in [11] and is illustrated in the next section in detail since it plays a crucial 
role in the understanding of the results of the simulation of cosmic rays in the Galaxy. 

The parallel diffusion coefficient can be estimated from the ID case, by using the 
quasi-linear theory: 

where J^{k) = {5B{k) / Bof is the normalized power in the turbulent modes with 
wavenumber /c oc 1/p resonant with the particles with momentum p. Even in the 
3D case this is a reasonable approximation to the parallel diffusion coefficient since this 
is dominated by the components of the perturbing field which are perpendicular to the 
background field. In this case one can see that the perpendicular diffusion coefficient is 

^ D\\J^{kf. (6) 

Since by definition T{k) <^ 1 it is easy to see that D_<^ <^ D\\, which implies that in 
most cases the effect of perpendicular diffusion is irrelevant if the propagation occurs in 
the regime of weak turbulence. 
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In numerical simulations of the propagation of cosmic rays in the Galaxy it is 
usually assumed that 6B/Bq ~ This ratio is supported by general estimates, 

such as equipartition, cosmic ray behavior and observations of total magnetic field in 
elliptical galaxies [16], rather than direct observations. 

Let us assume that the measurement of the abundances of light elements and the 
estimate of the anisotropy of cosmic rays at low energies may be taken as realistic for 
the determination of the diffusion properties of the ISM. 

The anisotropy of cosmic rays at low energies is observed to be at the level of 
S ~ 10"^, and in the context of quasi-linear theory (QLT) it is of order Vd/c, where Ve> 
is the drift velocity of cosmic rays in the magnetic field of the Galaxy. The condition 
vd/c ^ 10~^ implies ~ 3 x 10^ cm/s. This is in good agreement with the theory 
again, because in QLT the streaming instability forces the streaming of cosmic rays to 
occur at bulk velocities lower than the Alfven speed, va = B / y/^np ~ 2 x 10^ cm/s for 
B = and gas density O.lcm"^ (this should be considered as an average value over 
the magnetized halo of the Galaxy, say within 3 kpc from the disk) . In other words, the 
anisotropy is exactly what one would expect on the basis of bulk motion of cosmic rays 
at the Alfven speed {vd = va)- In QLT the pathlength for a particle to suffer a change 
in direction by 90 degrees is 



where Q = c/riij)) is the gyration frequency of the particle and k = l/r^^p). 
The pathlength A determines the confinement time in a region of size L as 

From observations of the abundance of light elements this time is measured to be 
~ 3 X 10^ years, while from the abundance of unstable radioactive isotopes one gets 
a larger number, ~ 2 x 10^ years ^7\. These two numbers correspond respectively to 
A = 10 pc and A = 1.5 pc. Here we assumed that the magnetized region of the Galaxy 
in the direction perpendicular to the disk has a typical size L = 3 kpc. Note also that 
rigorously we may use the parallel diffusion coefficient to estimate the escape from the 
disk only if the magnetic field is oriented along z, which is at odds with the conventional 
models of Galactic magnetic field. Therefore it is worth keeping in mind that a more 
realistic estimate is likely to differ from the one just illustrated and often used in the 
literature. 

From the equation for A one immediately obtains: 



A 




c 



(7) 




(9) 



for A = 10 pc and 




(10) 
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for A = 1.5 pc. For the numerical evaluation we considered cosmic rays with mean 
energy 1 GeV. These values of kP{k) correspond to 6B/Bq ~ 2 x 10~^ and ~ 5 x 10~^ 
respectively on the relevant scales. On such scales it appears that the assumptions of 
QLT are fulfilled. 

If the power spectrum is in the form of a power law, we can write P{k) = Pq (^^^ 
and limit ourselves to the two interesting cases a = 5/3 (Kolmogorov spectrum) and 
a = 3/2 (Kraichnan spectrum). In both these cases most power is in the form of 
modes with the largest spatial scale (namely at ko, assumed here to be ko ~ 1/100 pc). 
The modes of wavenumber ko resonate with particles with energy Eq = 2.8 x 10^'' 
eV. The propagation of particles with energies larger than Eq is described in terms of 
a diffusion coefficient with a steeper energy dependence than the one discussed here 
(Bohm diffusion) and eventually straight line propagation. From the numerical values 
obtained above, and assuming that ko ~ 1/100 pc, one easily infers that the power on a 
scale ko is 

Poko ~ ei (^^^^^)" ' = 3.2 X 10-^ (1.8 x 10"^) (11) 
for A = 10 pc and a = 5/3 {a = 3/2), and 

Poko ^ £2 (^^^^)" ' = 0-02 (1.3 X 10-3) (12) 

for A = 1.5 pc and a = 5/3 (a = 3/2). 

These estimates show that the total power in the turbulent field may be appreciably 
smaller than unity, which of course affects the normalization of the diffusion coefficient, 
the confinement time and the expected anisotropy at higher energies. The main problem 
with these estimates is that they are based solely upon the parallel diffusion coefficient, 
which, as discussed below may be incorrect. The issue of the strength of the turbulent 
field relative to the strength of the regular field remains therefore open. 

It is worth stressing that for a = 5/3 the diffusion approximation is broken at 
Eth ^ 8 X 10^^ eV when A = lOpc and Eth ^ 2 x 10^^ eV when A = 1.5 pc. For 
a = 3/2 we have Eth ~ 10^^ eV for A = 10 pc and Eth ~ 6 x 10^^ eV for A = 1.5 pc. 
This implies that at energy Eth the anisotropy is expected to become of order unity. 
Among all cases considered, the only case that seems to be compatible with the fact 
that no large anisotropy is observed up to the knee is the case a = 5/3 and A = 1.5 pc. 
Note that this does not necessarily imply that a large anisotropy should be observed 
at Eth ~ 2 X 10^^ eV, since at this energy the chemical composition in the Galaxy is 
expected to be contaminated by heavy elements, which are as isotropic as the particles 
with energy Eth/Z. Despite the interesting conclusion, this has to be considered just as 
a hint, because of the several assumptions that enter the previous estimate (for instance 
the value of L and ko and assumptions about geometry of the system). 

The predicted escape time from the Galaxy as a function of energy is more solidly 
predicted to be t{E) oc E~^^^ for Kolmogorov spectrum and t{E) oc E~^^'^ for Kraichnan 
spectrum. It is worth stressing that this simple prediction, widely used in the literature. 
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completely neglects the possibility of perpendicular diffusion or when it is not neglected, 
the assumption is adopted that the scaling with energy of the perpendicular diffusion 
coefficient is the same as for the parallel diffusion coefficient. Unfortunately, as we show 
below, the role of perpendicular diffusion in the Galaxy is likely to be crucial. 



3.2. Drifts 

Gradients in the modulus or orientation of the large scale field Bq induce drift motions 
in the direction perpendicular to both the local field and its gradient. The drift velocity 
of the guiding center can be written as [18]: 



cp r 1 . 2 -Bo X V5o , 2 ^0 X [{Bo ■ V)Bo] 

z^tr^" " Bi + " M 

1 2 ^0 X VBq 2 r^o X VBq (V X Bo 

cri\ - sm a h cos a 



Bi fin 



f J- . 2 -"0 ^ V ^(j 2 

where a is the pitch angle of the particle. The first term in this expression reflects the 
transverse gradient of the field strength while the second term represents the effect of 
the curvature of the field lines. 

The above expression should be interpreted as the drift velocity averaged over a 
gyration period of the particle. As an estimate of the order of magnitude of the time scale 
for escape from the region of size L due to drift motion, we can write T£){E) ~ crL{E) ' 
where \grad is the spatial scale on which the gradient in the magnetic field appears. This 
expression clearly shows that if the drifts are relevant at all this may happen only at 
very high energies. 

Three toy models are particularly interesting as far as drifts are concerned and will 
be discussed in detail in §5.21 §5.31 and §5.41 Here we limit our discussion to the expected 
effects of drift motions. The first model (Toy model II in §5.31) has only spatially constant 
(in modulus) azimuthal magnetic field. In this case the field lines are simply concentric 
circles in ^ = const, planes. The only gradient is due to the curvature of the magnetic 
field lines and the drift velocity is given by 

vr, = Eisc cos^ a —, (13) 
P 

where z is the unit vector in the z direction, p is the distance (in kpc) from the center in 
the plane z = and Eig is the particle energy in units of 10^^ eV. Clearly this expression 
and the ones we will list below are valid as long as the spatial scale of the gradient is 
much larger than the Larmor radius of the particles. This condition also assures that 
the drift velocity is always smaller than the speed of light. From Eq. [13] one can see 
that the drift pushes the particles perpendicular to the plane. 

The second toy model that we will consider is similar to the previous one but with 
the strength of the magnetic field having a gradient along the p direction (see Eq. 15. 3p . 
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It is easy to predict that also in this case the drift velocity is oriented along the z 
direction. The drift velocity in this case is 

vd = Eisc^{l + cos^ a) z for p > 4kpc. (14) 

Finally, in the third toy model we assume that the magnetic field is still azimuthal but 
is constant in the z = plane and has a gradient in the z direction (see Eq. I5.4p . In 
this case the drift velocity is 

V£, = Eisc cos^ a -i + sin^ a — p exp{z/zc). (15) 

L p 2Zc ^ 

Clearly in this third case the direction of the drift is no longer along z and depends on 
P- 

4. Determination of the parallel and perpendicular diffusion coefficients 

To calculate the diffusion coefficients we inject a few thousand particles of a given energy 
isotropically in a magnetic field composed of a constant regular component along z plus a 
uniform turbulent component. We record the particle trajectories and we then calculate 
the instantaneous parallel and perpendicular diffusion coefficients as: 

We plot the instantaneous diffusion coefficients as a function of propagation time in 
Fig. [T] for the case 6B/Bq = 1. The left panel is the parallel diffusion coefficient while the 
right one is the perpendicular diffusion coefficient. The parallel instantaneous diffusion 
coefficient increases linearly in the beginning when the particles are still only feeling 
the regular field and at some point flattens when the full diffusive regime is reached, 
typically within a few scattering lengths. In the perpendicular case the instantaneous 
diffusion coefficient increases for a time t-i/2 ^ tt x tl/c, corresponding to half a gyration 
around the regular field. At this point continuing the gyration the particle is going back 
to its starting position and the diffusion coefficient is decreasing, having a minimum at 
tl. After some time the diffusion regime is reached and the curve shows a plateau. This 
plateau identifies the diffusion coefficient and we use the average of the last 15 points 
(the gray region in the plots) to estimate it. 

In the following few paragraphs we present our results for the diffusion coefficients 
as a function of energy for some interesting configurations. 

4.1. The case of vanishing regular field 

Without a background field the only type of turbulence that can be considered is 
isotropic turbulence. In Fig. [2] we plot the diffusion coefficient as a function of energy 
for 3D turbulence in a configuration with no regular field, but only turbulent field with 
Linax = 100 pc and 6B = lOOpdJI. 

I Please note that here and in the foUowing when denoting SB = 100 fiG we actually mean: 
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Figure 1. Instantaneous diffusion cocfRcients as a function of propagation time. Left: 
parallel, right: perpendicular. Each line/color corresponds to a different particle energy 
as indicated in the plot. The black points on the far right of each line indicate the 
average of the corresponding points in the gray region and they represent the estimate 
of the diffusion coefficient at the corresponding energy. Note the different scales on 
y-axes. 



In this case, in order to compare our results with the ones of Ref. [H], we calculated 
the diffusion coefficients using 6 in the denominator of Eq. (fT6|) instead of 2. In Fig. [2] 
the gray points and lines are the diffusion coefficients along the three axes, the red 
points are the total diffusion coefficient and the black line is the parametrization of the 
diffusion coefficient given in Ref. [H] that was obtained from simulations using the plane 
wave approach. In this case we used the FFT approach and the agreement is very good. 

4.2. Combination of regular and turbulent fields 

In this case we use a superposition of a constant background field and a turbulent field 
with three levels of isotropic turbulence: 6B/Bq = 0.5,1,2. The maximum scale of 
the turbulence is set to Lmax = 0.1 kpc, Bq = 1 /iG and we use the FFT approach to 
generate the turbulence. We plot the parallel and perpendicular diffusion coefficients in 
Fig. [3l The top three lines represent the parallel diffusion coefficients, while the bottom 
three the perpendicular ones. The turbulence level is given by the numbers attached to 
the curves. It is interesting to note that while the low energy (10^^- 10^^ eV) slope of the 
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E[eV] 

Figure 2. Diffusion coefficient for a configuration with vanishing regular field. The 
gray points and lines are the diffusion coefficients along the three axes, the red points 
are the total diffusion coefficient and the black line is the parametrization of the 
diffusion coefficient given in Ref. [14] . 

parallel diffusion coefficient is 1/3 as one would expect, the slope of the perpendicular 
one is steeper, being about 0.5 — 0.6. 

In Fig. m we plot the ratio of the perpendicular to parallel diffusion coefficients, 
D±/D\\, as a function of energy. The three sets of points connected by solid lines 
are the results of the three simulations shown in the previous plot. We compared our 
results with the ones obtained in Ref. [H]. The thin black lines are the results from 
their Fig. 6 for the cases r] = 0.46 and 0.21 that correspond to SB/Bq = 0.92 and 0.52 
respectively. The agreement between the two sets of results is pretty good, especially 
for the case 5B/Bq ^ 1. The ratio of the diffusion coefficients is almost constant with 
a slow energy dependence. 

The tiny difference in slope between and D\\ at low energy is more apparent 
in this plot. It is interesting to note that this difference seems to be present also in 
the results of Ref. [H], at least for the case 6B/Bq ~ 1. In the case 6B/Bq ^ 0.5 the 
scattering of their points is too big to allow for inferring any conclusions in this respect. 
In order to confirm that this slope was not a systematic effect due to the method used 
to generate the turbulence, we performed a simulation using the plane wave approach 
to generate the field. These results are represented by the brown dashed line. The 
simulation parameters and the shape of the turbulence spectrum in this case are a bit 
different from the others, and the resulting curve does not coincide with the one obtained 
from the FFT approach, but also in this case the ratio is not constant at low energy 
and presents a small positive slope. 



Numerical propagation of high energy cosmic rays in the Galaxy I: technical issues 13 




^Q^5 ^Q^6 ^Q^7 

E[eV] 

Figure 3. Parallel and perpendicular diffusion coefficients as a function of energy for 
three levels of turbulence. The upper three lines are the parallel diffusion coefficients, 
while the bottom three represent the perpendicular one. The level of turbulence, 
SB/ Bq is given by the numbers attached to the lines. 

The results of Ref. |12] show no dependence of the ratio Dj^/ D\\ on energy and for 
6B/B0 = 1 their result is smaller than ours by about a factor 2. 

5. Toy models of the Galactic magnetic field 

The large scale structure of the Galactic magnetic field is likely to be complex, as 
made of spiral arms and various types of gradients along the radial direction in the 
disk and along the z axis, perpendicular to the disk. The same presence of the spiral 
arms induces gradients on different spatial scales. On top of this large scale structure 
a turbulent component is present which turns out to be responsible for the diffusive 
motion of cosmic rays. In all cases presented below, the values of the quantity 6B/Bq is 
assumed to be spatially constant (in other words the turbulent field is a constant fraction 
of the large scale field). It appears rather unrealistic that the naive expectations based 
on quasi-linear theory may find an easy confirmation with this complex structure of the 
magnetic field and indeed we confirm that this is in general not the case. In order to 
understand the various reasons why the expectations of QLT may be not fulfilled, in 
the following we discuss in detail four toy models of the magnetic field of the Galaxy 
in both its regular (large scale) and turbulent components. The first model is that of 
a magnetized homogeneous sphere with only turbulent field. In this case QLT cannot 
even be applied because of the absence of a regular field which does not allow to develop 
a perturbative approach to particle propagation. In this case however the confinement 




Figure 4. Ratio of tlie perpendicular to parallel diffusion coefficients, D±/D^^, as a 
function of the energy. The three set of points connected by solid lines are the results 
for the three levels of turbulence indicated. The points connected by the dashed line 
are the result of a simulation with a set of parameters similar to the one used for the 
orange one, but in this case using the plane waves approach instead of the FFT one. 
The two black thin lines are the results of Ref. [ij for SB/Bq = 0.92 (upper one) and 
0.52 (lower one). 



time that is obtained from simulations is close to the naive extrapolation of QLT to a 
regime in which it should not be applied. 

The second toy model consists of a purely azimuthal, spatially constant magnetic 
field. The particles are injected at the position of the Earth and collected on the surface 
of a cylinder of radius 10 kpc and height 0.5 kpc. 

The third and fourth toy models are variations of the previous one with the addition 
of gradients along the radial direction and along the z direction. 



5.1. Toy model I: a magnetized homogeneous sphere 

We consider a sphere filled uniformly with isotropic turbulent field with Lmax = 0.1 kpc 
and 6B = 0.5, l,2yuG. We inject protons in the center of the sphere and we collect 
them when they reach a distance of 2 kpc from the center. The times of escape from 
the sphere are plotted as triangles and boxes in Fig. [51 We also plotted, with stars, the 
results obtained using Lmax = 1 kpc instead of Lmax = 0.1 kpc for the case 6B = 1 fiG. 
The black lines are the expected propagation times obtained using the parametrization 
of the diffusion coefficient given in [T3] and already used in §4.11 for comparison: 
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Figure 5. Times of escape from a spiiere filled with uniform tm-bulcnt field for protons 
injected in the center. The levels of turbulence are indicated in the plot. The first 
three cases are for Lmax = 0.1 kpc, while the last one for i,nax = Ikpc. The black 
lines are the expected results obtained using the diffusion coefficient given in Ref . [Mj . 



The agreement is very good both in the low energy and in the transition region. Going 
to very high energies, the transition to straight hne propagation becomes visible. 

5.2. Toy model II: large scale azimuthal field with no spatial gradients 

The magnetic field as seen from above the disk is as shown in Fig. [6l This field structure 
is assumed to resemble at least qualitatively the spiral structure of the Galactic field. 
In passing we notice that this purely azimuthal field has also recently been adopted 
by [19]. The turbulent field is assumed to have a Kolmogorov spectrum with a largest 
scale Lmax = 0.1 kpc and total power 6B/Bq = 0.5,1 and 2. Particles are injected 
at the Earth, located at Rq = 8.5 kpc from the center and propagated backwards in 
time until they escape the cylinder of radius R = 10 kpc and height above and below 
the disk of 0.5 kpc. A crucial point to realize here is that the magnetic field lines are 
closed loops: the magnetic field strength is spatially constant but the orientation of the 
field changes as illustrated in Fig. [61 The fact that the field lines are closed implies 
a straightforward but important conclusion: the particles cannot escape the cylinder 
by diffusing parallel to the magnetic field lines. The only way particles can escape is 
by diffusing and drifting perpendicular to the field lines, which is clearly made more 
difficult by the smallness of the perpendicular diffusion coefficient (see §4j) as compared 
with the parallel diffusion coefficient. The escape times of cosmic rays as functions of 
energy for the various cases that have been calculated are illustrated in Fig. [7] (top 
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Figure 6. Azimuthal magnetic field in Toy model II. 

panel). The lower panel of the figure illustrates the column densities experienced by 
cosmic rays with given energy. The gas density has been assumed to be constant and 
equal to 1 cm~^ inside the disk (|2;| < 200 pc) and 0.01 cm~^ outside the disk. The 
different symbols refer to the values of 6B/Bq as indicated in the figure. The straight 
line represents the drift time calculated from Eq. ( fT3l) using the average drift velocity. 
It is worth noticing that the actual drift times have a very extended tail towards large 
times, due to the dependence of this quantity on the angle of injection of the particles 
with respect to the large scale local field. The black lines and dots are the diffusion 
time scales, oc D±{E)~^, where the perpendicular diffusion coefficient is taken from the 
simulations described in ^ 

A general comment about the relative role of parallel and perpendicular diffusion 
is in order: parallel diffusion is more effective than perpendicular even in the case of 
strong turbulence considered here, but it only leads to motion of the particles along the 
closed magnetic field lines. Perpendicular diffusion, though much slower, is responsible 
for particle escape in the direction perpendicular to the disk (there is also some escape 
from the sides of the cylinder but this process is less efficient because the sides are ~ 1.5 
kpc away from the location of the Earth, while the halo has been assumed to be only 0.5 
kpc thick). The parallel with the Galaxy is very instructive in this instance: particles 
diffuse effectively along the spiral arms, whose length is roughly Rarm. ~ t^Rq ~ 30 kpc 
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Figure 7. Particle escape times for Toy Model II. The upper panel shows the times 
required for the particles to escape from a cylinder with half-height of 0.5 kpc. The 
lower panel shows the grammage of gas crossed. The boxes and triangles are the values 
for different levels of turbulence as indicated in the plot. The injection is set at 8.5 kpc 
for all cases except for the light blue upward triangles for which it is set at 85 kpc. The 
thick black line is the drift timescale while the thin black curves represent the diffusion 
timescales. 
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long (at the distance of the Sun from the galactic center). The diffusion time parallel 
to the arms is therefore ry ~ R'^^^/D\\. At the same time, cosmic rays diffuse in the 
direction perpendicular to the disk in a time t± ~ R^/D±. The ratio of the two time 
scales is t\\/t± ~ 10^ D±/D^^, where we assumed ~ 1 kpc. For 6B/Bq ~ 1 the 
perpendicular diffusion coefficient is not much smaller than D\\, so that it is easy to 
understand that perpendicular diffusion may become the dominant channel of cosmic 
ray escape from the Galaxy rather than parallel diffusion. In our toy model this situation 
is extreme in that the magnetic field lines are closed and no escape at all is possible 
along the field. As a consequence, the energy dependence which is illustrated in Fig. [7] 
reflects the energy dependence of the perpendicular diffusion coefficient, which in the 
relevant energy range can be approximated as D± oc £'0.5-o.6 jg instructive that 
such a slope, usually associated (at low energies) to a Kraichnan spectrum of turbulence 
(parallel diffusion) can in fact be achieved with a Kolmogorov perpendicular diffusion 
(at least in the high energy range we are able to treat here). 

The important role of perpendicular diffusion in determining the escape time is also 
shown by the absolute normalization of the curves in Fig. [3 For parallel diffusion, at 
least in the quasi-linear regime, one expects the diffusion coefficient to decrease while 
increasing 6B/Bq, so that the escape times increase. In our toy model the perpendicular 
diffusion coefficient in fact increases with increasing 6B/Bq. 

Aside from diffusion, the escape times are also affected by drift motions. In 
particular, drifts become important where the drift time (the straight line in Fig. 
[7]) becomes of the same order of magnitude of the diffusion times (black lines). For 
SB/B = 0.5 this happens at energies around 10^^ eV, while drift seems irrelevant for 
stronger levels of turbulence. Besides this effect, which is rather clear from Fig. U\ 
there is a more subtle effect induced by drifts, which is evident in the low energy part 
of the curve for 6B/Bq = 0.5. One can notice that the black line illustrating the effect 
of diffusion (for SB/Bq = 0.5) lies below the upward red triangles obtained in the 
simulation. In order to understand the reason for this apparent problem, we calculated 
the escape times in the case in which the Earth is located at 85 kpc from the center 
instead of 8.5 kpc. In this case the cylinder is larger but it has the same height. But more 
important the curvature of the magnetic field lines is reduced appreciably so that the 
drift velocity drops correspondingly. One can see that the low energy behaviour in this 
case (upward light blue triangles) agrees well with the black curve, therefore confirming 
that the reason for the slim disagreement has something to do with the curvature of the 
field lines. 

To achieve a better understanding of the modifications the drifts produce to the 
diffusion process we calculated the diffusion coefficients in this geometry and we found 
that the drifts are modifying the two perpendicular diffusion coefficients reducing the 
one along z and increasing the one along p. In fact one should keep in mind that the 
concepts of parallel and perpendicular diffusion were introduced here with reference to 
the specific case of a large scale coherent background field, with no intrinsic curvature 
of the field lines. When the field lines are curved, then the definition itself of parallel 
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and perpendicular diffusion changes, as discussed in detail in Appendix A 

At energy ~ 10^^ eV the Larmor radius of the particles equals the maximum 
wavelength in the power spectrum of the turbulent field and the diffusion regime changes, 
gradually shifting toward the straight line propagation, which in Fig. [7| corresponds to 
the extreme right, flat part of the curves for the escape time. It is worth reminding the 
reader that all these considerations remain valid for heavier nuclei once the energy is 
substituted by rigidity. 

We conclude this discussion of the second toy model with a comment on the absolute 
magnitude of the escape time. Though keeping in mind that this is a toy model of the 
magnetic field of the Galaxy, we believe that some qualitative conclusions can be drawn. 
At energy 10^^ eV the escape time for the cases considered here is ris ^ 0.5 — 5 million 
years (the halo height here is only 0.5 kpc). These numbers are of the same order of 
magnitude of the confinement times estimated from the abundance of light element in 
the GeV region, which means that in order to fit these observations one should postulate 
that the escape time below 10^^ eV should be practically energy independent. We could 
not envision any realistic mechanism able to justify such an expectation. It follows 
that within the limitations of the toy model 2 it is very hard to obtain a realistic, even 
qualitative, description of what is observed in the Galaxy at much lower energies. This 
conclusion is confirmed also by the curves on the grammage: at 10^^ eV cosmic rays 
traverse a column density of 1 — 2 g cm~^. As we discuss below, this conclusion is the 
same for the other toy models considered here. 

5. 3. Toy model III: large scale azimuthal field with spatial gradient along p 

The global structure of the large scale azimuthal field is not changed with respect to 
Toy Model II, but we introduce here a gradient of the modulus of the large scale field 
with the radial coordinate p measured in the x — y plane. The radial dependence of the 
field is assumed to be in the form: 

' 2.125;uG p < 4 
^ 8.5/1^ p-^ p>4: ' 

where p is the radius in cylindrical coordinates in units of kpc. As discussed in §3.21 in 
this case the drift velocity is still oriented in the z direction, therefore the drift due to 
a gradient of the strength of the field behaves qualitatively as the gradient due to the 
curvature in the field lines, discussed in the section above. The escape time and the 
grammage for this case are illustrated in Fig. [HI where the red dashed line indicates the 
drift timescale, again calculated using the average drift velocity. At the distance of the 
Earth the gradient due to the radial dependence reduces the drift time by roughly a 
factor 2, thereby making the line for the drift time scale almost touch the red triangles 
{6B/Bq = 0.5). For stronger levels of turbulence the drifts become basically irrelevant, 
even at the highest energies. The absolute normalizations of the time scales are affected 
very little by the radial gradient of B, therefore most comments made for Toy Model 2 
are valid here too. 
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Figure 8. Particle escape times for Toy Model III. The black line is the drift timescale 
for the field of Toy Model II, while the red dashed line is the drift timescale for the 
present configuration. 



5.4- Toy model IV: large scale azimuthal field with spatial gradient along z 

We conclude this section by investigating the case of an azimuthal field with a gradient 
along £, as described by the following expression: 

B{z) = exp{—z/zc)fiG, 

with Zc = 0.25 kpc or 0.1 kpc. In this case the drifts due to the ^-dependence are in 
the radial direction and, at the Earth position, are bigger than the drifts due to the 
curvature of the field lines. The sum of the two drifts tends to push the particles toward 
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Figure 9. Particle escape times for Toy Model IV and Zc = 0.25 kpc. 



the center of the Galaxy, where the drifts due to curvature dominate. The exit points 
of the particles in this case are shifted in the direction of the galactic center, while in 
the previous two Toy Models most particles escaped from a ring with p ~ 8.5 kpc. 

The escape times and grammage for Zc = 0.25 kpc are illustrated in Fig. [9] with the 
usual meaning of the symbols. 

The effect of drifts, combined with the smaller effective size of the magnetized halo 
along z, contribute to reduce the escape times. At 10^^ eV the escape time is always 
shorter than 1 million year. However the slopes of the curves, although rather uncertain, 
do not seem to point toward any flattening that may help reconcile the grammage at 
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Figure 10. Particle escape times for Toy Model IV and Zc = 0.1 kpc. 



10^^ eV (0.2 — 1 g cm~^) with that observed in the GeV region. A further reduction of 
the escape times is achieved by reducing the scale Zc- For instance the time scales and 
grammage for Zc = 0.1 kpc are illustrated in Fig. [10 



6. Discussion and Conclusions 

The propagation of cosmic rays in the Galaxy still presents us with numerous open 
questions. The standard lore goes as follows: if the sources of galactic cosmic rays 
(possibly but not necessarily supernova remnants) inject a spectrum Q{E) oc E~"' with 
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7 ^ 2.1 — 2.4, then diffusion of these cosmic rays in the magnetic field of the Galaxy 
leads to an equilibrium spectrum which is n{E) oc E~"'~^, where the diffusion coefficient 
is taken in the form D{E) oc E^. For a Kolmogorov spectrum of magnetic fluctuations 
S = 1/3, while for a Kraichnan spectrum 6 = 1/2. Of course these statements apply at 
energies lower than the maximum energy of the accelerated particles, which for protons 
is expected to be ~ 10^^ — 10^^ eV. However, if the maximum energy of the accelerated 
particles were much larger, in principle the same conclusions would extend up to the 
energy for which the Larmor radius equals the coherence scale of the field, which is 
typically taken to be ~ 100 pc. This corresponds to energy ~ (1 — 3) x 10^'^ eV for 
a magnetic field 1 — 3fiG. The simulations illustrated in this paper can be performed 
for proton energies E > 10^^ eV (in a few cases E > 10^^ eV), therefore for at least 
two decades in energy we should be able to test the standard lore sketched above. We 
confirm that this is the case by considering a toy model with only a turbulent field with 
given power spectrum, in which case we are in perfect agreement with the expectations. 
The problems arise as soon as any complication is added to this simple scenario. We 
illustrate our points by considering other three toy models, each having a specific feature 
which is supposed to be resemblant of a corresponding feature expected to be present 
in the actual Galactic magnetic field. In particular we consider a benchmark situation 
in which the Galactic field is taken to have a perfectly azimuthal geometry, so that the 
magnetic field lines are closed loops. We showed how in such a geometry the role of 
perpendicular diffusion in the escape of particles from the toy Galaxy is crucial and 
leads to escape times which are too long to be reconciled with the observed confinement 
times at much lower energies. This conclusion should remain valid in the case in which 
the magnetic field lines follow the spiral arms rather than being closed, since the arms 
are in any case much larger in length than the size of the halo. 

An important piece of information should be added: the escape times that we 
plotted throughout the paper are all meant to be the average of the log of the escape 
times. The spread around these mean values are very large, covering about one order of 
magnitude. Such spreads do not reflect limitations in the statistics of the propagated 
particles: they are rather stable if the number of particles is increased. The fluctuations 
are due to the several possible histories that may characterize the propagation of cosmic 
rays in the Galaxy. On the other hand, the mean values used to infer our conclusions 
are very stable. 

Another important ingredient of the magnetic field configuration with closed 
magnetic field lines consists of the drift motions induced on the particles by the gradients 
in the direction of the local large scale field. The effect of drift is especially evident for 
high energies and low levels of turbulence. Similar drifts are introduced by gradients in 
the z and p directions. 

The most important conclusion that we could achieve is that the dominant role of 
the perpendicular diffusion in a geometry with a prominent azimuthal magnetic field 
makes the expectation of the common lore hard to realize. The energy dependence of 
the perpendicular diffusion coefficient is not the same as that of the parallel diffusion 
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coefficient: more specifically in the lower energy regime it scales as D± oc ^E'O-s-o-s^ 
rather than E^^^ as would be expected for a Kolmogorov spectrum. Unfortunately we 
are not able to follow this behaviour down to energies below 10^^ eV. In any case, at 
~ 10^^ eV, the escape times that we measured in the simulation are always too large 
to be extrapolated down to the few million years inferred from the abundance of light 
elements in the GeV energy region, even admitting that a flattening to a behaviour 
oc E~^^'^ of the escape times could be achieved below ~ 10^^ eV. 

It is interesting to speculate about possible physical effects that might cause the 
escape time to be reduced. From the discussion above, it is clear that reducing the 
level of turbulence (namely the value of 6B/B) does not help, since this would cause 
the perpendicular diffusion to decrease, thereby increasing the escape times even more. 
Making the halo have a smaller scale height does help, but it appears rather unrealistic to 
reduce this scale below 0.5 kpc (observations of the radio background from synchrotron 
emission of relativistic electrons hint to a typical scale height of a few kpc [20]). One 
possibility that we will discuss more quantitatively in Paper II is that of a galactic wind, 
possibly injected by cosmic rays themselves: in this case, in addition to the diffusive 
motion, particles would have a systematic drift velocity pushing them away from the 
disc of the Galaxy. If the typical wind velocity is uw ~ 10 ''cm s~^ [21], then the typical 
escape time due to advection is of order 5 million years, independent of momentum. It 
is important to notice that for the diffusion coefficients used in the literature (in the 
common lore usually one does not distinguish between parallel and perpendicular) this 
is roughly the escape time scale for cosmic rays in the GeV energy region, therefore the 
effect of the wind is usually relevant only for low energy particles (at higher energies 
the escape is dominated by diffusion). In the scenarios that we find, escape is due to 
perpendicular diffusion, and takes place on much longer time scales as we have seen, 
therefore the effect of the wind can be that of producing a roughly energy independent 
escape time of the order of ~ 5 million years. Unfortunately this does not appear to be 
the correct, or at least the complete, picture either. In fact the escape time is observed to 
be a function of energy r oc E~^-^, as shown by the energy dependence of the secondary 
to primary ratio (e.g. [22]), although these measurements have so far been carried out 
only up to energies of the order of 10^ — 10^ MeV/nucleon. 
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Appendix A. Calculation of the diffusion coefficients in the azimuthal field 

We consider a regular magnetic field with constant magnitude and azimuthal direction 
as in §5.21 We inject particles at po = 8.5 kpc or po = 85kpc and we propagate them 
for 1 Mpc recording their trajectories. In this case we cannot use Eq. [16] to calculate the 
diffusion coefficients since the average values of the displacements we are considering 
are no longer due to the drifts. 

We proceed as follows: we build histograms of the particle positions at fixed times 
during the propagation and then we fit these histograms with gaussian distributions. 
The fitted value of the variance allows us to estimate the diffusion coefficient while the 
mean value of the distribution is related to the drift velocity. 

The three directions we used to calculate the diffusion coefficients are: z, p and 
(f). The first two correspond to the two perpendicular coefficients and the latter to the 
parallel one. For z we simply histogram the z coordinate of the particles. For p we 
histogram the p = y/x'^ + y'^ and then we divide each bin in the histogram by p to take 
into account the volume element (in cylindrical coordinates). For we histogram x po, 
where po is the distance at which the particles were injected. 

In Fig. lAll we plot the three diffusion coefficients as a function of propagation 
time for injection at 8.5 kpc and 6B/B = 0.5. The top panel represents the diffusion 
coefficient in the p direction, the middle panel the one in the z direction and the lower 
panel the parallel diffusion coefficient. The differently colored lines represent, from 
bottom to top, different energies from 10^^ eV to 10^^'^ eV with a logarithmic step of 
0.2. The points to the far right of the plots are the average of the last 10 corresponding 
points and represent our estimate of the diffusion coefficient. 

Concerning the parallel diffusion coefficient we can see that the curves are flat 
and that the diffusion regime is achieved. The only unexpected feature is in the two 
highest energy curves, corresponding to 10^^ and 10^^'^ eV that show a steepening around 
cr ~ 1 Mpc. This steepening is simply due to the fact that at high energies and large 
propagation times some of the particles have enough time to complete half a circle 
around the "galaxy" and since to measure the parallel displacement we are using x po , 
particles with |0| > vr end up in the wrong place in the histogram and distort the 
distribution. This is not however a physical effect and it is just a glitch of the method 
used to estimate the parallel displacement and we can just trow away the last few points 
and do the average with the remaining ones. 

The diffusion coefficient along z shows a tiny sub-diffusion at low energies, 
Dz{E, t) oc r~°'^^, that disappears by increasing the particle energy. It is interesting to 
note that in the plots of Fig. [U that were obtained for similar parameters, but with the 
large scale field constant and along the z direction, the diffusion regime was obtained 
already with cr ~ 10 kpc, with slightly larger times necessary for higher energies. In 
the present situation the results show that the opposite is occurring: at high energy 
the particles reach the diffusion regime, whereas at low energy they may not, at given 
time. In this case, since the instantaneous diffusion coefficient shows sub-diffusion, it is 
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Figure Al. Instantaneous diffusion coefficients as a function of propagation time for 
injection at 8.5 kpc in a field composed of large scale azimuthal field and an isotropic 
turbulent field with SB/Bq — 0.5. The three panels, from top to bottom, show the 
three diffusion coefficients along p, z and (f> respectively. 



not completely correct to define a diffusion coefficient using the average of the last few 
points. We do it anyway averaging the points with propagation times between 100 kpc 
and 1 Mpc that represent the range of propagation times obtained in Fig. [7| for energies 
between 10^^ eV and 10^ ''eV. In this way we obtain at least a rough estimate of the 
diffusion coefficient affecting the particle propagation in our specific case. 

For the diffusion coefficient in the p direction we have a situation similar to the 
z one, but with super-diffusion instead of sub-diffusion. In this case the effect is even 
smaller with: Dp{E,T) oc r°'^. Again increasing the energy the anomalous diffusion 
disappears. 
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Increasing the injection distance to 85 kpc the anomalous diffusion is reduced, but it 
is still slightly present. On the other hand increasing the turbulence level to SB/Bq = 1 
or 2 reduces it much more than increasing the distance. 

This is clear in the plots of Fig. IA2I where we report the diffusion coefficients as 
a function of energy for the three levels of turbulence and the two injection distances. 
The black thick lines are the results for the case of constant large scale field directed 
along z (the curves of Fig. [3]). The green lines are the diffusion coefficients in the (j) 
direction, the parallel ones. The red and blue lines are the diffusion coefficients in the p 
and z direction respectively. The solid lines are for injection at 8.5 kpc and the dotted 
ones for injection at 85 kpc. 

The above caveat about anomalous diffusion notwithstanding, the plots in Fig. IA2I 
seem to explain the results of Fig.[7l For the case of injection at 8.5 kpc and SB/Bq = 0.5 
we found in §5.2l that the obtained escape time was bigger that the one we expected from 
the perpendicular diffusion coefficient. This is consistent with the results presented in 
Fig. IA2I where it is shown that, in this case, the diffusion coefficient in the z direction 
is reduced and this obviously produces an increase in the escape time. Increasing the 
injection distance the z diffusion coefficient is closer to the "unmodified" one (see dotted 
blue line in the top panel of Fig. IA2p and the times of escape are almost on top of the 
expectations (see light blue triangles and top black line in Fig. [7]). 

Increasing the level of turbulence, the effect of the curvature of the field lines is 
reduced and both the p and z diffusion coefficients rapidly converge towards the normal 
one (middle and bottom panels in Fig. IA2[) . 
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Figure A2. Diffusion coefficients as a function of energy. Green lines: parallel. Red 
lines: along p. Blue lines: along z. Black lines: parallel and perpendicular diffusion 
coefficients from Fig. [S] Solid lines: injection at 8.5 kpc. Dotted lines: injection at 
85 kpc. The three levels of turbulence used are indicated in the panels. 
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